close all;
clear all;
l0=0.1; % rest length m
m = 0.1; % 1 kg mass
k = 1500; %spring stiffness
b = 50;%m/seconds
tMax = 1; %seconds
deltaT = 0.001; %seconds
Fapp = 25; %Newtons

x1Min =l0/4;
x2Min=2*x1Min; 

x1 = l0 ; % initial extention
xDot1 =0; %initial velocity
xDoubleDot1 = 0;

x2= 2*l0 ; % initial extention
xDot2 =0; %initial velocity
xDoubleDot2 = 0;

t= 0;
store1=x1;
store2=x2; 

Fapp1 = 0;
Fapp2 = 0;
%Fapp2 = Fapp;
%Fapp = 0;

while (t<tMax)

	%disp(t);	
	if (( t > 0.10000) && (Fapp1!=Fapp))
		Fapp1=Fapp;
		disp("Fapp1 has changed");
	endif
	
	xDoubleDot1 = (1/m)*( - Fapp1 - ( k*(x1-l0)) - (b*xDot1) + k*(x2-(x1+l0)) );

	%Integration of Positionn and velocity;
	xDot1 = xDot1 + xDoubleDot1*deltaT;
	x1 = x1 + xDot1*deltaT;
%	if (x1<x1Min)
%		x1=x1Min;
%	endif
		
	% Changing the Force on mass 2 after half the time		
	if ((t>(0.5)*(tMax)) && (Fapp2!=Fapp))
		disp("Changing Fapp2");
		Fapp2 = Fapp;
	endif
	
	xDoubleDot2 = (1/m)*(- Fapp2 - (k * (x2-(l0+x1))) - (b*xDot2));

	%Integration of Position and Velocity
	xDot2 = xDot2 + xDoubleDot2*deltaT;
	x2 = x2 + xDot2*deltaT;
	
%	if (x2<x2Min)
%		x2=x2Min;
%	endif

	t= t+deltaT;
	store1= [store1;x1];
	store2= [store2;x2];
	
	plot(x1, 5, 'o-', x2, 5, 'og'); 
	axis([0,0,10,10],"square");
	pause(0.001);
endwhile

%plot (store1, 'r');  hold on; 
%plot(store2 , 'b'); 
